/* 
 * File:   LU.h
 * Author: gustavo
 *
 * Created on June 17, 2012, 12:04 PM
 */

#ifndef LU_HPP
#define	LU_HPP

#include "../matrix.hpp"

#include "LUInterface.hpp"

namespace imp {

	namespace math {

		template <typename T, template <typename> class Algorithm > class LU : public LUInterface<T> {
		private:
			int pivsign;

		public:
			Matrix<T> M;

		public: /* Constructors */
			LU(const Matrix<T>& orig);
			virtual ~LU();

		public: /* Implements LUInterface<T> */
			void getLower(Matrix<T>& Lower);
			void getDiagonal(Matrix<T>& Diagonal);
			void getUpper(Matrix<T>& Upper);
			void getPermutation(Matrix<T>& Permutation);
			bool isSingular();
			T determinant();
			Matrix<T> solve(Matrix<T>& B);

		public: /* Implements DecompositionInterface< Matrix<T> > */
			bool modifiedInput();
		private:
			void decompose(Matrix<T>& orig);

		};

		/* = Constructors ======================================================= */
		template <typename T, template <typename> class Algorithm > LU<T, Algorithm>::LU(const Matrix<T>& orig) {
			M.copy(orig);
			this->decompose(M);
		}

		template <typename T, template <typename> class Algorithm > LU<T, Algorithm>::~LU() {

		}

		/* = Implements LUDecomp<T> ============================================= */
		template <typename T, template <typename> class Algorithm > void LU<T, Algorithm>::getLower(Matrix<T>& Lower) {

			int rows = M.rows;
			int cols = M.rows < M.columns ? M.rows : M.columns;

			Lower.setSize(rows, cols);
			//Lower.clear();

			for (int i = 0; i < cols; i++) {
				Lower(i, i) = 1;

				for (int j = 0; j < i; j++) {
					Lower(i, j) = M(i, j);
				}
			}

			if (rows > cols) {
				for (int i = cols; i < rows; i++) {
					for (int j = 0; j < cols; j++) {
						Lower(i, j) = M(i, j);
					}
				}
			}
		}

		template <typename T, template <typename> class Algorithm > void LU<T, Algorithm>::getDiagonal(Matrix<T>& Diagonal) {

		}

		template <typename T, template <typename> class Algorithm > void LU<T, Algorithm>::getUpper(Matrix<T>& Upper) {

		}

		template <typename T, template <typename> class Algorithm > void LU<T, Algorithm>::getPermutation(Matrix<T>& Permutation) {

		}

		template <typename T, template <typename> class Algorithm > bool LU<T, Algorithm>::isSingular() {
			for (int i = 0; i < M.rows; i++) {
				if (fabs(M(i, i)) < imp::internal::EPS) {
					return true;
				}
			}
			return false;
		}

		template <typename T, template <typename> class Algorithm > T LU<T, Algorithm>::determinant() {

			if (M.rows != M.columns) {
				__error("Must be a square matrix.\n");
				//throw new Exception("Must be a square matrix.");
			}

			double ret = pivsign;

			for (int i = 0; i < M.size; i += M.columns + 1) {
				ret *= M.data[i];
			}

			return ret;
		}

		template <typename T, template <typename> class Algorithm > Matrix<T> LU<T, Algorithm>::solve(Matrix<T>& B) {

		}

	}

}

#endif	/* LU_HPP */

